################################################
# Clean replication code for Tandeo 2022: Descriptive plots
# Created AJH
# Created 5/4/2022
# Last edited 5/4/2022
################################################
library(tidyverse)
library(srvyr)
library(survey)
library(stargazer)
library(sf)
library(scales)
library(gridExtra)

# Plot of the survey boroughs
alc <- read_sf("chapter_5_rewarding_responsiveness/additional_data/alcaldias.shp")
alc$in_sample <- ifelse(alc$cve_mun %in% c("005", "006", "007", "011", "017"),1,0)
# note that I manually annotate this after the fact

p <-  ggplot(alc) + geom_sf(aes(fill = factor(in_sample)))+ 
  theme_bw()+
  theme(legend.position="none")+
  labs(title = "Survey Data Is From the City's\n Five Eastern Alcaldias")   +
  scale_fill_grey(start=.8, end = .2)
ggsave(filename = "/Users/alyssahuberts/Dropbox/2_mx_water/1_Tandeo/3_Output/1_Plots/5_Spring2022/graphics/survey_map.jpg",p)


  # load colonia-level tandeo data
  load("chapter_5_rewarding_responsiveness/additional_data/tandeo_colonias.Rdata")
  # use only 2009b
  tandeo_colonias$tandeo_2009a <- NULL
  tandeo_colonias <- tandeo_colonias %>%  rename(tandeo_2009=tandeo_2009b)
  rtypes <- tandeo_colonias %>%
    pivot_longer(cols =c(tandeo_2004, tandeo_2005, tandeo_2006, tandeo_2007, tandeo_2008, tandeo_2009, tandeo_2010, 
                       tandeo_2011, tandeo_2012, tandeo_2013, tandeo_2014, tandeo_2015, tandeo_2016, tandeo_2017, 
                       tandeo_2018, tandeo_2019, tandeo_2020, tandeo_2021), 
               names_to= "year", 
               values_to ="tandeo",
               names_prefix="tandeo_") %>% 
  group_by(cve_col) %>% summarise(years_in_sample=n(),
                                                  years_rationed = sum(tandeo),
                                                  percent_rationed = years_rationed/years_in_sample)

  colonias_sf <- read_sf("/Users/alyssahuberts/Dropbox/2_mx_water/1_Tandeo/2_Data/9_Administrative/coloniascdmx/coloniascdmx.shp")
  colonias_sf <- left_join(colonias_sf, rtypes, by="cve_col" )

  p1 <-  ggplot(colonias_sf)+
      geom_sf(aes(fill=percent_rationed), lwd=.05)+
      theme_classic()+
      scale_fill_distiller(type = "seq",
                            direction = 1,
                            palette = "Greys") +
      theme(legend.position = "bottom")+
      labs(title = "Historical Inclusion on\n Rationing Lists",
           subtitle = "Data from 2004 to 2021",
           x= "",
           y = "",
           fill = "Percent of Period \n Rationed")+
    geom_sf_text(data = colonias_sf[colonias_sf$cve_col=="15-037",], label = "Historic \n Center", size = 2)

      ggsave(plot = p1, file = "/Users/alyssahuberts/Dropbox/2_mx_water/1_Tandeo/3_Output/1_Plots/5_Spring2022/graphics/years_rationed_continuous.jpg")

    # read in altitude data
      alt <- read_csv("chapter_5_rewarding_responsiveness/additional_data/altitude_seccion.csv")
      alt$cve_secc <- str_pad(alt$SECCION, 4, "left","0")
      secciones_sf <- read_sf("chapter_5_rewarding_responsiveness/additional_data/secciones.shp") %>% select(cve_secc, geometry)
      alt_sf <- left_join(secciones_sf, alt, by = "cve_secc")
      p2 <- ggplot(alt_sf)+
        geom_sf(aes(fill=alt_mean), lwd=.05)+
        theme_classic()+
        scale_fill_distiller(type = "seq",
                             direction = 1,
                             palette = "Greys", breaks= c(seq(1900, 3700,300)), labels=comma) +
        theme(legend.position = "bottom")+
        labs(title = "Altitude",
             x= "",
             y = "",
             fill = "Average Altitude (KM) \n Precinct")

        

      
      ggsave(plot = p2, file = "/Users/alyssahuberts/Dropbox/2_mx_water/1_Tandeo/3_Output/1_Plots/5_Spring2022/graphics/altitude.jpg")
     
      p <- grid.arrange(p2,p1, ncol=2) 
      ggsave(plot = p, file = "/Users/alyssahuberts/Dropbox/2_mx_water/1_Tandeo/3_Output/1_Plots/5_Spring2022/graphics/altitude_and_rationing.jpg")
      
      
      # Calculate census percentage with access to piped water
      census <- read_csv("chapter_2_making_scarcity_work/principales_resultados_localidad.csv") %>% clean_names
      census <- census %>%  filter(nom_mun == "Total de la entidad Ciudad de México"& nom_loc == "Total de la Entidad")
      as.numeric(census$vph_aguadv)/as.numeric(census$tvivparhab)

